function DecodingProbability_v27(q,g,T,save_var)

% Clear commandline
clc
close all

% Time it
tic

% Transmitter range
T_min=T(1);
T_step=T(2);
T_max=T(3);

% Field cardinality is the q vector input
% Layer selection probability is the g vector input

% Layer dimensions (these are dejan examples)
k1=32;
k2=64;
K1=k1;
K2=k1+k2;

% Create all transition matrices
global PP_q1 PP_q2
PP_q1=cell(1,K2);
PP_q2=cell(1,K2);
for i=K1:K2
    PP_q1{i}=TransitionMatrix(i,q(1));
    PP_q2{i}=TransitionMatrix(i,q(2));
end

disp('Done with all transition matrices')

% Calculate decoding probabilities #set 1
l1_prob_1=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q(1),g);
l2_prob_1=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q(1),g);

% Calculate decoding probabilities #set 2
if save_var==1
    l1_prob_2=layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q(2),g);
    l2_prob_2=layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q(2),g);
end

% Save data

% Save data, plot
if save_var == 1
    
    % Plot nice graph!
    hold all
    
    % Linespec options here
    % http://www.mathworks.se/help/techdoc/ref/linespec.html
    
    plotter(l1_prob_1,l2_prob_1,T_min,T_step,T_max,K1,K2,q(1));
    legend_h=plotter(l1_prob_2,l2_prob_2,T_min,T_step,T_max,K1,K2,q(2));
    
    % Set proper legend
    label_1=strcat('L1, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^1)');
    label_2=strcat('L2, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^1)');
    label_3=strcat('L1, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^8)');
    label_4=strcat('L2, \Gamma=',num2str(g(1)),' , ',num2str(g(2)),', FF(2^8)');
    set(legend_h,'String',[label_1;label_2;label_3;label_4],'location','NorthWest')
    set(legend_h,'interpreter','tex')
    
    figname=strcat('uep_ew_analytic','_g1_0',num2str(g(1)*10),'_g2_0',num2str(g(2)*10),'.eps');
    print(gcf,'-deps',figname)
    workspace=strcat('uep_ew_analytic','_g1_0',num2str(g(1)*10),'_g2_0',num2str(g(2)*10));
    save(workspace);
else
    % Write values to terminal
    disp('The probability of decoding layer 1 with q=2 is:')
    l1_prob_1(T_min)
    save('EWcomparisondata')
end



% Time it
toc

end

% Calculate layer 1 probability
function l1_prob = layer_1_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

l1_prob=zeros(T_max,1); % Decoding 1. Layer probabilities

for tx =T_min:T_step:T_max % For each number of recv packets
    
    sol=0;
    
    for n = 0:tx % For each permutation of a number of recv packets
        
        % Layer 1 by itself
        val=PM2(n,K1,K1,q);
        
        val_test=0;
        
        % Layer 1 and Layer 2 gives rank K2 (This way we also get L1!)
        % Sum prob for all possible ways to achieve rank K2 with given permutation of recv packets
        
        for i=0:K1-1
             
            tmp1=PM2(n,K1,i,q);   
            
            %Optimization We only need to calculate tmp2 if tmp1!=0
            if tmp1==0
                continue
            end
            
            tmp2=PM2(tx-n,K2-i,K2-i,q);
            val_test=val_test+tmp1*tmp2;
            
        end
        
        % We have counted all the ways L1 can become full rank by itself
        % We have counted all the ways L2 can become full rank (except when
        % L1 is full rank!)
        % Both outcomes are valid for getting L1 and since they are disjoint we can just add them!
        % P(a)+P(b)-P(ab), where P(ab)=0 because they are disjoint!
        
        sol=sol+(val+val_test)*binopdf(n,tx,g(1));
        
    end
    
    l1_prob(tx)=sol;
    disp(['Layer 1: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Calculate layer 2 probability
function l2_prob = layer_2_decod_prob(T_min,T_step,T_max,K1,K2,q,g)

% Decoding 2. Layer probabilities
l2_prob=zeros(T_max,1);

for tx =T_min:T_step:T_max % For each number of recv packets
    
    sol=0;
    
    for n = 0:tx % For all permutations of recv packets
        
        val=0;
        
        for i=0:K1 % For all possible ways to achieve rank K2 with given permutation of recv packets
            
            tmp1=PM2(n,K1,i,q);
            
            % Optimization no need for tmp2 when tmp1=0!
            if tmp1==0
                continue;
            end
            
            tmp2=PM2(tx-n,K2-i,K2-i,q);
            val=val+tmp1*tmp2;
        end
        
        sol=sol+val*binopdf(n,tx,g(1));
        
    end
    
    l2_prob(tx)=sol;
    
    disp(['Layer 2: ' num2str(tx) ' out of ' num2str(T_max)])
    
end

end

% Plotter for a nice graph!
function legend_h = plotter(l1_prob,l2_prob,T_min,T_step,T_max,K1,K2,q)

% Replace 0 with NaN in (l1_prob,l2_prob) for prettier plot
for k=1:length(l1_prob)
    
    if l1_prob(k)==0
        l1_prob(k)=NaN;
    end
end

for k=K2:length(l2_prob)
    if l2_prob(k)==0
        l2_prob(k)=NaN;
    end
    
end


% Plotting
figure(1)

if q ==2^1
    l1_style='s';
    l2_style='--s';
else
    l1_style='o';
    l2_style='--o';
end

% Legend fix
plot(-1,-1,'-s',-1,-1,'--s',-1,-1,'-o',-1,-1,'--o');

plot(1:length(l1_prob),l1_prob,'-','Color','k','MarkerSize',4)
plot(4:4:length(l1_prob),l1_prob(4:4:end),l1_style,'Color','k','MarkerSize',4)

plot(1:length(l2_prob),l2_prob,'--','Color','k','MarkerSize',4)
plot(4:4:length(l2_prob),l2_prob(4:4:end),l2_style,'Color','k','MarkerSize',4)

% Plot annotation
grid('on')
pbaspect([2.5 1 1])
legend_h = legend('location','NorthWest');
xlabel('Total number of received packets [-]')
ylabel('Decoding probability [-]')
set(gca,'XTick',0:10:T_max)
set(gca,'YTick',0:0.1:1)
xlim([T_min T_max])
ylim([0 1])

end

% The New helper function
function P = PM2(m,n,r,q)

global PP_q1 PP_q2

s1=zeros(n+1,1);
s1(1)=1;

if q==2^1
    val=(PP_q1{n}^m)*s1;
    P=val(r+1);
end

if q==2^8
    val=(PP_q2{n}^m)*s1;    
    P=val(r+1);
end


end

% Returns a transition matrix of dim 'n+1' by 'n+1' with q param
function M = TransitionMatrix(n,q)
P1=zeros(n+1,n+1);
for i=1:length(P1)
    P1(i,i)=1/(q^(n-(i-1)));
    if i<n+1
        P1(i+1,i)=1-P1(i,i);
    end 
end

M=P1;

end















